###   Implementation Constraints on natural climate solutions
###   Author: Kroeger, Timm; Erbaugh, James; Luo, Zhixian; Hilary Brumberg; Waverly Eichhorst; Margaret Hegwood; Anna LoPresti; Priya Shyamsundar; Peter W. Ellis; Lauren E. Oakes; Dow Martin; Pedro H. S. Brancalion; Mieke Bourne; Arundhati Jagadish; Kemen G. Austin; Andrew Kinzer; Marcos Sanjuán; Lisa McCullough; Marta Echavarria,

library(tidyverse)
library(ggplot2)
library(grid)

nbase<-read.csv("NBaseData_v1.csv")

### Figure 1. Constraint Total
constrainttotals <- nbase %>%
  group_by(Constraint_category, Source_type) %>%
  summarize(count = n()) %>%
  ungroup() %>%
  complete(Constraint_category, Source_type, fill = list(count = 0)) %>%
  group_by(Constraint_category) %>%
  mutate(total_count = sum(count)) %>%
  ungroup()

## Generate the % lable for each category
labels_figure1_df <- constrainttotals %>%
  dplyr::distinct(Constraint_category, total_count) %>%
  dplyr::mutate(pct_lab = (total_count / sum(total_count))) %>%
  dplyr::mutate(pct_lab = round(pct_lab * 100, 0)) %>%
  dplyr::mutate(pct_lab = paste0(pct_lab, "%"))
  


# Create the plot
totals_projcnstrnt <- ggplot(constrainttotals, aes(x = fct_reorder(Constraint_category, total_count, .fun = sum), y = count, fill = Constraint_category)) +
  geom_col(aes(fill = interaction(Constraint_category, Source_type)), position = "stack") +
  geom_text(data = labels_figure1_df,
            aes(x = fct_reorder(Constraint_category, total_count, .fun = sum),
                y = total_count, label = pct_lab),
            hjust = -0.08, size = 5 / .pt, color = "grey20", family = "Helvetica") +
  coord_flip(clip = "off") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
  theme_classic() +
  ylab(NULL) + xlab(NULL) + labs(fill = "") +
  theme(text = element_text(family = "Helvetica"),
        axis.text.x = element_text(size = 5, hjust = 1),
        axis.text.y = element_text(size = 5),
        legend.position = "none",
        axis.line  = element_line(linewidth = 0.06 * .pt, colour = "black"),
        axis.ticks = element_line(linewidth = 0.06 * .pt, colour = "black"),
        axis.ticks.length = unit(0.8, "pt")) +
  scale_fill_manual(values = c("Finance.S" = "#77aadd", 
                               "Finance.L" = "#77aadd80",
                               "Knowledge.S" = "#99ddff",
                               "Knowledge.L" = "#99ddff80",
                               "Material Inputs.S" = "#44bb99",
                               "Material Inputs.L" = "#44bb9980",
                               "Social and Behavioral.S" = "#bbcc33",
                               "Social and Behavioral.L" = "#bbcc3380",
                               "Government and Organizations.S" = "#aaaa00",
                               "Government and Organizations.L" = "#aaaa0080",
                               "Markets.S" = "#eedd88",
                               "Markets.L" = "#eedd8880",
                               "Negative side effects of NCS.S" = "#ee8866",
                               "Negative side effects of NCS.L" = "#ee886680",
                               "Rules and Laws.S" = "#ffaabb",
                               "Rules and Laws.L" = "#ffaabb80"))

totals_projcnstrnt

ggsave("Fig 1.pdf", totals_projcnstrnt,  width=88, height=60, units = "mm", dpi=300)

## Figure 2. Individual constraint totals
subconstrainttotals <- nbase %>%
  group_by(Constraint_category, Source_type, Constraint_short) %>%
  summarize(count = n()) %>%
  ungroup()

labels_figure2_df <- subconstrainttotals %>%
  group_by(Constraint_short) %>%
  summarize(total_count = sum(count)) %>%
  ungroup() %>%
  dplyr::mutate(pct_lab = (total_count / sum(total_count))) %>%
  dplyr::mutate(pct_lab = round(pct_lab * 100, 2)) %>%
  dplyr::mutate(pct_lab = paste0(pct_lab, "%"))


#####   Visualizing subconstraint frequencies   #####
totals_projsubcnstrnt <- ggplot(subconstrainttotals, aes(x = fct_reorder(Constraint_short, count, .fun = sum, .desc = FALSE), y = count, fill = interaction(Constraint_category, Source_type))) +
  geom_text(data = labels_figure2_df,
            aes(x = fct_reorder(Constraint_short, total_count, .fun = sum),
                y = total_count, label = pct_lab),
            hjust = -0.08, size = 5 / .pt, color = "grey20",
            inherit.aes = FALSE,
            family = "Helvetica") +
  geom_col(position = "stack") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
  coord_flip(clip = "off") +
  theme_classic() +
  ylab(NULL) + xlab(NULL) + labs(fill = "") +
  theme(text = element_text(family = "Helvetica"),
        axis.text.x = element_text(size = 5, hjust = 1),
        axis.text.y = element_text(size = 5),
        legend.position = "none",
        axis.line  = element_line(linewidth = 0.06 * .pt, colour = "black"),
        axis.ticks = element_line(linewidth = 0.06 * .pt, colour = "black"),
        axis.ticks.length = unit(0.8, "pt")) + # Remove the legend
  scale_fill_manual(values = c("Finance.S" = "#77aadd", 
                               "Finance.L" = "#77aadd80",
                               "Knowledge.S" = "#99ddff",
                               "Knowledge.L" = "#99ddff80",
                               "Material Inputs.S" = "#44bb99",
                               "Material Inputs.L" = "#44bb9980",
                               "Social and Behavioral.S" = "#bbcc33",
                               "Social and Behavioral.L" = "#bbcc3380",
                               "Government and Organizations.S" = "#aaaa00",
                               "Government and Organizations.L" = "#aaaa0080",
                               "Markets.S" = "#eedd88",
                               "Markets.L" = "#eedd8880",
                               "Negative side effects of NCS.S" = "#ee8866",
                               "Negative side effects of NCS.L" = "#ee886680",
                               "Rules and Laws.S" = "#ffaabb",
                               "Rules and Laws.L" = "#ffaabb80"))

totals_projsubcnstrnt

legend_labels <- c("Finance", "Knowledge", "Material Inputs", "Social and Behavioral", "Government and Organizations", "Markets", "Negative side effects of NCS", "Rules and Laws")
legend_colors <- c("#77aadd", "#99ddff", "#44bb99", "#bbcc33", "#aaaa00", "#eedd88", "#ee8866", "#ffaabb")


legend_grob <- grobTree(
  textGrob("Legend", x = unit(1, "npc") - unit(1, "cm"), y = unit(0, "npc") + unit(10, "cm"), hjust = 1, vjust = 0, gp = gpar(fontsize = 5, fontfamily = "Helvetica")), # Moved to the top
  rectGrob(x = unit(1, "npc") - unit(1, "cm"), y = unit(0, "npc") + unit(9, "cm"), width = unit(0.5, "cm"), height = unit(0.5, "cm"), gp = gpar(fill = legend_colors[1], col = NA)),
  textGrob(legend_labels[1], x = unit(1, "npc") - unit(1.5, "cm"), y = unit(0, "npc") + unit(9, "cm"), hjust = 1, vjust = 0, gp = gpar(fontsize = 5, fontfamily = "Helvetica")), # Adjusted x position
  rectGrob(x = unit(1, "npc") - unit(1, "cm"), y = unit(0, "npc") + unit(8, "cm"), width = unit(0.5, "cm"), height = unit(0.5, "cm"), gp = gpar(fill = legend_colors[2], col = NA)),
  textGrob(legend_labels[2], x = unit(1, "npc") - unit(1.5, "cm"), y = unit(0, "npc") + unit(8, "cm"), hjust = 1, vjust = 0, gp = gpar(fontsize = 5, fontfamily = "Helvetica")), # Adjusted x position
  rectGrob(x = unit(1, "npc") - unit(1, "cm"), y = unit(0, "npc") + unit(7, "cm"), width = unit(0.5, "cm"), height = unit(0.5, "cm"), gp = gpar(fill = legend_colors[3], col = NA)),
  textGrob(legend_labels[3], x = unit(1, "npc") - unit(1.5, "cm"), y = unit(0, "npc") + unit(7, "cm"), hjust = 1, vjust = 0, gp = gpar(fontsize = 5, fontfamily = "Helvetica")), # Adjusted x position
  rectGrob(x = unit(1, "npc") - unit(1, "cm"), y = unit(0, "npc") + unit(6, "cm"), width = unit(0.5, "cm"), height = unit(0.5, "cm"), gp = gpar(fill = legend_colors[4], col = NA)),
  textGrob(legend_labels[4], x = unit(1, "npc") - unit(1.5, "cm"), y = unit(0, "npc") + unit(6, "cm"), hjust = 1, vjust = 0, gp = gpar(fontsize = 5, fontfamily = "Helvetica")), # Adjusted x position
  rectGrob(x = unit(1, "npc") - unit(1, "cm"), y = unit(0, "npc") + unit(5, "cm"), width = unit(0.5, "cm"), height = unit(0.5, "cm"), gp = gpar(fill = legend_colors[5], col = NA)),
  textGrob(legend_labels[5], x = unit(1, "npc") - unit(1.5, "cm"), y = unit(0, "npc") + unit(5, "cm"), hjust = 1, vjust = 0, gp = gpar(fontsize = 5, fontfamily = "Helvetica")), # Adjusted x position
  rectGrob(x = unit(1, "npc") - unit(1, "cm"), y = unit(0, "npc") + unit(4, "cm"), width = unit(0.5, "cm"), height = unit(0.5, "cm"), gp = gpar(fill = legend_colors[6], col = NA)),
  textGrob(legend_labels[6], x = unit(1, "npc") - unit(1.5, "cm"), y = unit(0, "npc") + unit(4, "cm"), hjust = 1, vjust = 0, gp = gpar(fontsize = 5, fontfamily = "Helvetica")), # Adjusted x position
  rectGrob(x = unit(1, "npc") - unit(1, "cm"), y = unit(0, "npc") + unit(3, "cm"), width = unit(0.5, "cm"), height = unit(0.5, "cm"), gp = gpar(fill = legend_colors[7], col = NA)),
  textGrob(legend_labels[7], x = unit(1, "npc") - unit(1.5, "cm"), y = unit(0, "npc") + unit(3, "cm"), hjust = 1, vjust = 0, gp = gpar(fontsize = 5, fontfamily = "Helvetica")), # Adjusted x position
  rectGrob(x = unit(1, "npc") - unit(1, "cm"), y = unit(0, "npc") + unit(2, "cm"), width = unit(0.5, "cm"), height = unit(0.5, "cm"), gp = gpar(fill = legend_colors[8], col = NA)),
  textGrob(legend_labels[8], x = unit(1, "npc") - unit(1.5, "cm"), y = unit(0, "npc") + unit(2, "cm"), hjust = 1, vjust = 0, gp = gpar(fontsize = 5, fontfamily = "Helvetica")) # Adjusted x position
)

# Save the plot with the specified size and resolution
ggsave("Fig 2.pdf", plot = totals_projsubcnstrnt + annotation_custom(legend_grob), width = 180, height = 180, units = "mm", dpi = 300)

## Figure 3. Z-score heatmap by NCS pathway

nbase_figure4 <- nbase %>%
  mutate(NCS_Pathway = ifelse(NCS_Pathway == "regenerative agriculture (other than agroforestry)", "regenerative agriculturen\n (other than agroforestry)", NCS_Pathway)) %>%
  mutate(NCS_Pathway = ifelse(NCS_Pathway == "reduced woodfuel harvest in forests", "reduced woodfuel harvest\n in forests", NCS_Pathway)) %>%
  mutate(NCS_Pathway = ifelse(NCS_Pathway == "avoided coastal wetland conversion", "avoided coastal\n wetland conversion", NCS_Pathway)) %>%
  mutate(NCS_Pathway = ifelse(NCS_Pathway == "avoided grassland conversion", "avoided grassland\n conversion", NCS_Pathway))

nbase_heatmap <- nbase_figure4 %>%
  group_by(Constraint_category, NCS_Pathway) %>%
  summarise(frequency = n()) %>%
  ungroup()

heatmapscaled <- nbase_heatmap %>%
  group_by(NCS_Pathway)%>%
  summarise(avg = mean(frequency), stDev = sd(frequency)) %>%
  full_join(nbase_heatmap,.)%>%
  mutate(., zscore=((frequency-avg) / stDev))

allcombinations <- expand.grid(NCS_Pathway = unique(heatmapscaled$NCS_Pathway), Constraint_category = unique(heatmapscaled$Constraint_category))
heatmapscaled <- merge(heatmapscaled, allcombinations, by = c("NCS_Pathway", "Constraint_category"), all.y = TRUE)

scaled_heatmap<-ggplot(heatmapscaled, aes(x=NCS_Pathway, y=Constraint_category, fill=zscore))+
  geom_tile() +
  scale_fill_gradient(name="Z-Scores", low = "white", high = "red", na.value = "grey")+
  coord_flip()+
  theme_classic()+
  ylab(NULL)+xlab(NULL)+
  theme(text = element_text(family = "Helvetica", size = 5),
        legend.text = element_text(size = 5),
        axis.text.x = element_text(size = 5, angle = 45, hjust = 1),
        axis.text.y = element_text(size = 5),
        axis.line  = element_line(linewidth = 0.06 * .pt, colour = "black"),
        axis.ticks = element_line(linewidth = 0.06 * .pt, colour = "black"),
        axis.ticks.length = unit(0.8, "pt"))

scaled_heatmap

ggsave("Fig 3.pdf", scaled_heatmap,  width=88, height=88, units = "mm", dpi = 300)

## Figure 4. Z-score heatmap by subregion
nbase_heatmap <- nbase %>%
  group_by(Constraint_category, Subregion) %>%
  summarise(frequency = n()) %>%
  ungroup()

heatmapscaled <- nbase_heatmap %>%
  group_by(Subregion)%>%
  summarise(avg = mean(frequency), stDev = sd(frequency)) %>%
  full_join(nbase_heatmap,.)%>%
  mutate(., zscore=((frequency-avg) / stDev))

allcombinations <- expand.grid(Subregion = unique(heatmapscaled$Subregion), Constraint_category = unique(heatmapscaled$Constraint_category))
heatmapscaled <- merge(heatmapscaled, allcombinations, by = c("Subregion", "Constraint_category"), all.y = TRUE)


scaled_heatmap<-ggplot(heatmapscaled, aes(x=Subregion, y=Constraint_category, fill=zscore))+
  geom_tile() +
  scale_fill_gradient(name="Z-Scores", low = "white", high = "red", na.value = "grey")+
  coord_flip()+
  theme_classic()+
  ylab(NULL)+xlab(NULL)+
  theme(text = element_text(family = "Helvetica", size = 5),
        legend.text = element_text(size = 5),
        axis.text.x = element_text(size = 5, angle = 45, hjust = 1),
        axis.text.y = element_text(size = 5),
        axis.line  = element_line(linewidth = 0.06 * .pt, colour = "black"),
        axis.ticks = element_line(linewidth = 0.06 * .pt, colour = "black"),
        axis.ticks.length = unit(0.8, "pt"))
scaled_heatmap

ggsave("Fig 4.pdf", scaled_heatmap,  width=88, height=88, units = "mm", dpi = 300)